HDT2
Marisa Montoya, Majo Morales, Luis Garcia y Juanfer De Leon
8/2/2021
##Analisis exploratorio
Consumo Dataset
## 'data.frame': 257 obs. of 27 variables:
## $ X : int 3 4 5 6 7 8 9 10 11 12 ...
## $ Anio : int 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
## $ Mes : int 1 2 3 4 5 6 7 8 9 10 ...
## $ GasolinaSuper : num 308157 307766 331910 315648 319668 ...
## $ GasolinaRegular : num 202645 205531 229500 210680 208164 ...
## $ TotalGasolinas : num 510802 513297 561410 526328 527832 ...
## $ Diesel : num 634667 642381 699807 586804 656948 ...
## $ DieselLS : num NA NA NA NA NA NA NA NA NA NA ...
## $ DieselULS : num NA NA NA NA NA NA NA NA NA NA ...
## $ TotalDiesel : num 634667 642381 699807 586804 656948 ...
## $ GLP : num 194410 174711 189234 174331 191745 ...
## $ GasolinaAviacion : num 1426 1458 1503 1561 1642 ...
## $ Kerosina : num 64026 62660 61362 61814 54098 ...
## $ TurboJet : num NA NA NA NA NA NA NA NA NA NA ...
## $ Bunker : num 296767 328116 368590 396300 449369 ...
## $ Asfalto : num 48446 50597 27593 53794 60137 ...
## $ PetCoke : num NA NA NA NA NA NA NA NA NA NA ...
## $ AceitesLubricantes: num NA NA NA NA NA NA NA NA NA NA ...
## $ GrasasLubricantes : num NA NA NA NA NA NA NA NA NA NA ...
## $ Solventes : num NA NA NA NA NA NA NA NA NA NA ...
## $ Naftas : num NA NA NA NA NA NA NA NA NA NA ...
## $ Ceras : num NA NA NA NA NA NA NA NA NA NA ...
## $ CrudoNacional : num NA NA NA NA NA NA NA NA NA NA ...
## $ Butano : num NA NA NA NA NA NA NA NA NA NA ...
## $ Orimulsion : int NA NA NA NA NA NA NA NA NA NA ...
## $ MezclasOleosas : num NA NA NA NA NA NA NA NA NA NA ...
## $ Total : num 1750545 1773220 1909499 1800933 1941772 ...Podemos ver que todas las variables tratadas en las datasets son variables cuantitativas, por lo que haremos histogramas de las variables.
## consumo$Anio :
## Frequency Percent Cum. percent
## 2020 12 4.7 4.7
## 2019 12 4.7 9.3
## 2018 12 4.7 14.0
## 2017 12 4.7 18.7
## 2016 12 4.7 23.3
## 2015 12 4.7 28.0
## 2014 12 4.7 32.7
## 2013 12 4.7 37.4
## 2012 12 4.7 42.0
## 2011 12 4.7 46.7
## 2010 12 4.7 51.4
## 2009 12 4.7 56.0
## 2008 12 4.7 60.7
## 2007 12 4.7 65.4
## 2006 12 4.7 70.0
## 2005 12 4.7 74.7
## 2004 12 4.7 79.4
## 2003 12 4.7 84.0
## 2002 12 4.7 88.7
## 2001 12 4.7 93.4
## 2000 12 4.7 98.1
## 2021 5 1.9 100.0
## Total 257 100.0 100.0
Variable X: no sigue distribución normal, y es lo esperado ya que X es una variable diferente para los datos, como un ID.
Variable Anio: No sigue una distribucion normal tampoco, pero nos muestra que hay una frecuencia mayor en los primeros años comenzando en 2000, luego las frecuencias en los demas años son similares hasta 2021, cuando los datos son muy pocos.
Variable mes: No sigue una distribución normal tampoco.
Variable Gasolina Super: No sigue una distribución normal.
Variable Gasolina Regular: No sigue una distribución normal.
Variable Gasolina Diesel: Si sigue una distribución normal, por lo que pueda que tenga un comportamiento peculiar en el estudio de las series lineales.
####Cruce de variables Al estar trabajando con las variables de gasolina, en los cruces de variables los haremos siempre con las variables de gasolinaSuper, Regular y Diesel.
De este cruce de variables podemos ver que puede que exista un tipo de tendencia, que con los años, el consumo de gasolina super creció con los años.
Podemos ver que lo mismo sucedió con gasolina Regular y los años, pero otra vez, es la diesel la que tiene un comportamiento irregular.
Ahora vamos a hacer el cruce de variables pero con meses.
plot(y=consumo$GasolinaSuper,x=consumo$Mes,ylab="Gasolina Super", xlab="Mes")
plot(y=consumo$GasolinaRegular,x=consumo$Mes,ylab="Gasolina Regular", xlab="Mes")
plot(y=consumo$TotalDiesel,x=consumo$Mes,ylab="Gasolina Diesel", xlab="Mes")
Ahora vemos que las variables se comportan de una manera aleatoria conforme a los meses del año. Ya sea cualquier tipo de las gasolinas tratadas.
Analisis Importacion
Breve resumen del data set de importaciones
## [1] "Dimensiones"## [1] 245 29## [1] "Resumen"## X Anio Mes GasolinaSuper
## Min. : 3.0 Min. :2001 Min. : 1.000 Min. : 170293
## 1st Qu.: 64.0 1st Qu.:2006 1st Qu.: 3.000 1st Qu.: 357350
## Median :261.0 Median :2011 Median : 6.000 Median : 435416
## Mean :267.5 Mean :2011 Mean : 6.429 Mean : 464691
## 3rd Qu.:433.0 3rd Qu.:2016 3rd Qu.: 9.000 3rd Qu.: 558403
## Max. :623.0 Max. :2021 Max. :12.000 Max. :1227174
##
## GasolinaRegular TotalGasolinas Diesel DieselLS
## Min. : 81015 Min. : 274500 Min. : 229765 Min. : 691066
## 1st Qu.:194830 1st Qu.: 558276 1st Qu.: 620918 1st Qu.: 891341
## Median :282509 Median : 718114 Median : 767285 Median :1056569
## Mean :341982 Mean : 806673 Mean : 782290 Mean :1056457
## 3rd Qu.:469437 3rd Qu.:1021878 3rd Qu.: 904976 3rd Qu.:1195728
## Max. :896841 Max. :1861582 Max. :1595699 Max. :1592580
## NA's :41 NA's :204
## DieselULS TotalDiesel GLP GasolinaAviacion
## Min. : 1203 Min. : 229765 Min. :100562 Min. : 48.0
## 1st Qu.: 7176 1st Qu.: 655715 1st Qu.:210871 1st Qu.: 487.2
## Median :17701 Median : 800223 Median :371098 Median : 2132.0
## Mean :21905 Mean : 829959 Mean :375817 Mean : 3292.6
## 3rd Qu.:38300 3rd Qu.: 994338 3rd Qu.:510507 3rd Qu.: 3367.5
## Max. :48946 Max. :1630636 Max. :960841 Max. :27979.1
## NA's :225 NA's :89
## Kerosina TurboJet Bunker Asfalto
## Min. : 1000 Min. : 18373 Min. : 8485 Min. : 74.88
## 1st Qu.: 33762 1st Qu.: 45061 1st Qu.: 172130 1st Qu.: 2996.07
## Median : 51860 Median : 69728 Median : 294609 Median : 5497.36
## Mean : 54898 Mean : 73330 Mean : 308321 Mean : 7231.09
## 3rd Qu.: 67657 3rd Qu.: 95020 3rd Qu.: 417863 3rd Qu.:10237.67
## Max. :184094 Max. :158719 Max. :1051764 Max. :48364.50
## NA's :50 NA's :183 NA's :10
## PetCoke AceitesLubricantes GrasasLubricantes Solventes
## Min. : 6993 Min. :15086 Min. :115.8 Min. : 3207
## 1st Qu.:149221 1st Qu.:18020 1st Qu.:313.8 1st Qu.: 6264
## Median :165003 Median :21592 Median :388.2 Median : 8881
## Mean :218887 Mean :22826 Mean :436.4 Mean : 9974
## 3rd Qu.:244990 3rd Qu.:24292 3rd Qu.:546.0 3rd Qu.:11978
## Max. :849463 Max. :44746 Max. :996.5 Max. :24501
## NA's :127 NA's :216 NA's :216 NA's :216
## Naftas Ceras Butano PetroleoReconstit
## Min. : 0.380 Min. :107.9 Min. : 0.37 Min. :356364
## 1st Qu.: 3.263 1st Qu.:297.5 1st Qu.: 72.29 1st Qu.:360619
## Median : 54.000 Median :404.4 Median : 72.29 Median :367573
## Mean :138.722 Mean :428.8 Mean : 76.75 Mean :489360
## 3rd Qu.:280.017 3rd Qu.:522.2 3rd Qu.: 81.30 3rd Qu.:718114
## Max. :350.920 Max. :920.1 Max. :149.09 Max. :730957
## NA's :233 NA's :218 NA's :206 NA's :225
## MTBE Orimulsion MezclasOleosas PetroleoCrudo
## Min. : 8184 Min. :249983 Min. : 166.7 Min. :191
## 1st Qu.:11062 1st Qu.:312030 1st Qu.: 966.5 1st Qu.:191
## Median :12680 Median :315916 Median :1668.6 Median :191
## Mean :12952 Mean :315469 Mean :1626.9 Mean :191
## 3rd Qu.:14090 3rd Qu.:325083 3rd Qu.:2417.1 3rd Qu.:191
## Max. :19431 Max. :344685 Max. :3134.8 Max. :191
## NA's :234 NA's :232 NA's :223 NA's :244
## TotalMensual
## Min. :1381787
## 1st Qu.:2065597
## Median :2468706
## Mean :2558891
## 3rd Qu.:2924855
## Max. :4322010
## ## [1] "Estructura"## 'data.frame': 245 obs. of 29 variables:
## $ X : int 3 4 5 6 7 8 9 10 11 12 ...
## $ Anio : int 2001 2001 2001 2001 2001 2001 2001 2001 2001 2001 ...
## $ Mes : int 1 2 3 4 5 6 7 8 9 10 ...
## $ GasolinaSuper : num 373964 243091 312084 285055 300914 ...
## $ GasolinaRegular : num 177777 123116 161726 127339 168730 ...
## $ TotalGasolinas : num 551740 366207 473811 412394 469644 ...
## $ Diesel : num 566102 489526 575560 437745 552609 ...
## $ DieselLS : num NA NA NA NA NA NA NA NA NA NA ...
## $ DieselULS : num NA NA NA NA NA NA NA NA NA NA ...
## $ TotalDiesel : num 566102 489526 575560 437745 552609 ...
## $ GLP : num 194066 170703 161837 163049 171519 ...
## $ GasolinaAviacion : num 820 3054 677 3399 585 ...
## $ Kerosina : num 33834 67440 31787 25801 45529 ...
## $ TurboJet : num NA NA NA NA NA NA NA NA NA NA ...
## $ Bunker : num 214582 294609 315264 205653 278371 ...
## $ Asfalto : num 27749 7504 26304 7886 8443 ...
## $ PetCoke : num NA NA NA NA NA NA NA NA NA NA ...
## $ AceitesLubricantes: num NA NA NA NA NA NA NA NA NA NA ...
## $ GrasasLubricantes : num NA NA NA NA NA NA NA NA NA NA ...
## $ Solventes : num NA NA NA NA NA NA NA NA NA NA ...
## $ Naftas : num NA NA NA NA NA NA NA NA NA NA ...
## $ Ceras : num NA NA NA NA NA NA NA NA NA NA ...
## $ Butano : num NA NA NA NA NA NA NA NA NA NA ...
## $ PetroleoReconstit : int 715344 370166 360530 359527 723346 360369 360570 719303 360635 717717 ...
## $ MTBE : int 8402 NA NA 8184 12680 NA 10642 13357 NA 19431 ...
## $ Orimulsion : num NA NA NA NA NA NA NA NA NA NA ...
## $ MezclasOleosas : num NA NA NA NA NA NA NA NA NA NA ...
## $ PetroleoCrudo : int NA NA NA NA NA NA NA NA NA NA ...
## $ TotalMensual : num 2312639 1769209 1945770 1623638 2262727 ...Aqui se encuentran los histogramas de las variables cuantitativas para ver si hay normalidad aun que sea en las continuas. 
Podemos observar que las variables: GasolinaSuper, TotalGasolinas, Diesel y TotalDiesel siguen una distribución normal.
Cruze de variables con el anio y mes 
En las tres gráficas anteriores del cruce de variables de las importaciones de cualquier tipo de gasolina y el año, se puede observar que hay una tendencia a crecer. Esto implica que conforme pasa el tiempo, cada vez se importa más gasolina sin importar su tipo.
Por otro lado, al observar las gráficas del cruce de variables de cualquier tipo de gasolina con los meses, no se puede observar alguna tendencia. Esto implica que sin importar el mes, siempre va a haber una importación de combustible más o menos constante.
Para comenzar con el analisis de series de tiempo nos encontramos con la necesidad de unir la columna de mes y anio en una sola columna en formato de fecha, luego crear diferentes data frames para los 3 tipos de gasolinas para hacer una linea del tiempo de cada una.
####Linea de tiempo para consumo de gasolina super Primero se convierte el data frame en una serie de tiempo. Luego podemos observar el inicio, fin, frecuencia y vista de la linea del tiempo.
## [1] 2000 1## [1] 2021 5## [1] 12
De esta linea podemos ver a primera vista que posee una tendencia a crecer, y que hay momentos en la linea de tiempo que hace que las varianzas tambien se vean afectadas y sean diferentes entre si. Podemos ver que se trata de una serie con frecuencia anual donde hay datos de todos los meses desde enero de 2000 hasta mayo del 2021.Como se puede ver la cantidad de galones vendidos va aumentando a medida que pasan los años.
#####Descomposicion de la serie 
Podemos observar una serie con tendencia a aumentar, que no es estacionaria en varianza, y además tiene estacionalidad.
Ahora vemos los datos de inicio y fin del train y del test
## [1] 2000 1## [1] 2014 12## [1] 2015 1## [1] 2021 5#####Estimar los parámetros del modelo. Cómo no es estacionaria en varianza le haremos una transformación logaritmica para hacerla constante en varianza.
Veremos si logramos estacionarizar la serie con esta transformación:
Podemos observar que las varianzas si se pudieron hacer mas constantes, sin embargo, se pueden observar picos de varianza desiguales, lo cual puede ser un problema y concreta que no se pudo estacionarizar en varianza en su totalidad.
Por lo que ahora trataremos de quitar la tendencia y estacionarizar en normalidad.
##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: 0.3673
## P VALUE:
## 0.7319
##
## Description:
## Fri Aug 06 21:42:55 2021 by user: jfdel##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: 0.3673
## P VALUE:
## t: 0.7895
## n: 0.773
##
## Description:
## Fri Aug 06 21:42:55 2021 by user: jfdelAmbas pruebas de Dickey-Fuller Test nos demuestra que p es mayor a 0.05 por lo que no podemos rechazar Ho, y por ende, asumir que no hay raices unitarias. Por lo que la serie de tiempo no es estacionaria en media.
Por lo que debemos de hacer lo mismo pero con una diferenciacion en la serie de tiempo.
## Warning in adfTest(diff(train)): p-value smaller than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -17.2389
## P VALUE:
## 0.01
##
## Description:
## Fri Aug 06 21:42:55 2021 by user: jfdelComo se puede ver el valor de p ahora sà está por debajo de 0.05 por lo que se puede descartar la hipótesis nula de que existan raices unitarias.
Solo se hizo uso de una diferenciación, por lo que d es 1, pero para conocer p y q se debe de realizar un gráfico de autocorrelación y autocorrelación parcial.
acf(logGS,50) 
pacf(logGS,50) 
Del ACF podemos ver que la gráfica no se anula y tiene un decrecimiento rapido. pero en el PACF se anula luego de 1, por lo que q es 1. P no se anula en ningun momento. por lo que proponemos un modelo, podemos pensar que p es 4.
Por lo que podemos proponer un modelo ARIMA con p =4 q=1 d = 1
Veamos si hay estacionalidad en la serie. Si volvemos a ver la descomposición de la serie en su componente estacional.
Al parecer sà existe estacionalidad en la serie. Para tener una idea de los parámetros estacionales veremos las funciones de autocorrelación y autocorrelación parcial con 36 resagos para ver en que momentos son los picos estacionales. Se usará la serie estacionarizada.
Cómo podemos ver en la función de autocorrelación tenemos un decaimiento en los picos estacionales, 12, 18, 24, 30, 36, y en la función de autocorrelación parcial, un pico significativo en el retardo 12. Eso sugiere que los parámetros del componente estacional son P = 2, D = 1 (el perÃodo es de 6 meses por lo que se puede hacer una diferenciación estacional), y Q=1.
Creamos los analisis del modelo Hacemos un analisis de coeficientes del AutoArima y Arima.
## Warning: package 'lmtest' was built under R version 4.0.5##
## Attaching package: 'lmtest'## The following object is masked from 'package:epiDisplay':
##
## lrtest##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.530999 0.255698 -2.0767 0.0378324 *
## ma1 -0.016694 0.240298 -0.0695 0.9446149
## ma2 -0.462657 0.139363 -3.3198 0.0009008 ***
## sar1 -0.021098 0.080054 -0.2635 0.7921275
## sar2 -0.137446 0.080886 -1.6993 0.0892688 .
## sma1 -0.999970 0.089828 -11.1321 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.694455 0.061480 -11.2957 < 2.2e-16 ***
## sma1 0.403711 0.079137 5.1014 3.371e-07 ***
## sma2 0.218524 0.067869 3.2198 0.001283 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vemos que todos los coeficientes son significativos y vamos a hacer un analisis de residuos.
Análisis de Residuales Los residuales deben ser lo más parecidos posible a un ruido blanco. Esto es que sean aleatorios (no están correlacionados entre sÃ), y que estén distribuidos normalmente. Analizamos los residuos de nuestro modelo fitArima
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,2)(2,1,1)[12]
## Q* = 13.882, df = 18, p-value = 0.7367
##
## Model df: 6. Total lags used: 24Según los gráficos podemos ver que la distribución de los datos parece ser normal, y que no hay correlaciones significativas. Según el test de Ljung-Box los datos se distribuyen de forma independiente puesto que el p-value es mayor a 0.05, por lo que no se puede rechazar la hipótesis nula. Esto significa que el modelo generado es aceptable para predecir.
Analizando los residuos del modelo generado de forma automática por R tenemos los siguientes resultados:
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,0,2)[12]
## Q* = 35.203, df = 21, p-value = 0.02682
##
## Model df: 3. Total lags used: 24El modelo autogenerado no es un modelo aceptable para predecir.
##Predicción con el modelo generado
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 13.03072 12.96888 13.09256 12.93614 13.12530
## Feb 2015 13.00463 12.93677 13.07249 12.90085 13.10841
## Mar 2015 13.09947 13.02944 13.16950 12.99237 13.20658
## Apr 2015 13.04365 12.96995 13.11735 12.93094 13.15636
## May 2015 13.02189 12.94554 13.09824 12.90512 13.13866
## Jun 2015 12.97006 12.89072 13.04940 12.84872 13.09139
## Jul 2015 13.04983 12.96784 13.13183 12.92443 13.17523
## Aug 2015 13.05805 12.97337 13.14273 12.92854 13.18755
## Sep 2015 13.01030 12.92308 13.09752 12.87690 13.14370
## Oct 2015 13.04912 12.95939 13.13885 12.91189 13.18635
## Nov 2015 13.02783 12.93569 13.11997 12.88691 13.16875
## Dec 2015 13.16829 13.07378 13.26281 13.02374 13.31284
## Jan 2016 13.04320 12.94569 13.14071 12.89407 13.19233
## Feb 2016 13.02564 12.92563 13.12564 12.87270 13.17858
## Mar 2016 13.11409 13.01176 13.21642 12.95758 13.27059
## Apr 2016 13.06203 12.95736 13.16669 12.90196 13.22209
## May 2016 13.03888 12.93197 13.14579 12.87537 13.20239
## Jun 2016 13.00052 12.89139 13.10965 12.83362 13.16743
## Jul 2016 13.07097 12.95968 13.18227 12.90076 13.24119
## Aug 2016 13.07571 12.96229 13.18914 12.90224 13.24919
## Sep 2016 13.02441 12.90889 13.13992 12.84774 13.20107
## Oct 2016 13.06679 12.94922 13.18437 12.88698 13.24661
## Nov 2016 13.04360 12.92401 13.16318 12.86071 13.22648
## Dec 2016 13.17061 13.04903 13.29219 12.98468 13.35655
## Jan 2017 13.05617 12.93330 13.17904 12.86826 13.24408
## Feb 2017 13.03266 12.90822 13.15709 12.84235 13.22296
## Mar 2017 13.12479 12.99868 13.25089 12.93193 13.31764
## Apr 2017 13.08114 12.95346 13.20883 12.88586 13.27642
## May 2017 13.05784 12.92856 13.18712 12.86012 13.25556
## Jun 2017 13.01275 12.88191 13.14360 12.81265 13.21286
## Jul 2017 13.08970 12.95731 13.22210 12.88723 13.29218
## Aug 2017 13.09458 12.96065 13.22850 12.88976 13.29939
## Sep 2017 13.03977 12.90434 13.17521 12.83264 13.24690
## Oct 2017 13.09205 12.95512 13.22898 12.88263 13.30147
## Nov 2017 13.06462 12.92622 13.20303 12.85295 13.27630
## Dec 2017 13.20159 13.06171 13.34147 12.98766 13.41552
## Jan 2018 13.08123 12.93921 13.22324 12.86404 13.29842
## Feb 2018 13.05667 12.91300 13.20033 12.83695 13.27638
## Mar 2018 13.14960 13.00440 13.29479 12.92754 13.37166
## Apr 2018 13.10526 12.95849 13.25203 12.88080 13.32972
## May 2018 13.08215 12.93386 13.23045 12.85536 13.30895
## Jun 2018 13.03535 12.88554 13.18517 12.80623 13.26448
## Jul 2018 13.11345 12.96213 13.26477 12.88202 13.34487
## Aug 2018 13.11880 12.96598 13.27161 12.88509 13.35250
## Sep 2018 13.06455 12.91027 13.21884 12.82860 13.30051
## Oct 2018 13.11613 12.96038 13.27188 12.87794 13.35433
## Nov 2018 13.08906 12.93187 13.24625 12.84866 13.32946
## Dec 2018 13.22767 13.06903 13.38630 12.98505 13.47028
## Jan 2019 13.10596 12.94511 13.26681 12.85996 13.35196
## Feb 2019 13.08224 12.91973 13.24475 12.83371 13.33077
## Mar 2019 13.17465 13.01061 13.33869 12.92378 13.42552
## Apr 2019 13.12917 12.96356 13.29478 12.87589 13.38245
## May 2019 13.10608 12.93894 13.27322 12.85046 13.36170
## Jun 2019 13.06024 12.89157 13.22891 12.80228 13.31820
## Jul 2019 13.13742 12.96724 13.30760 12.87716 13.39768
## Aug 2019 13.14274 12.97106 13.31442 12.88018 13.40530
## Sep 2019 13.08896 12.91581 13.26212 12.82414 13.35379
## Oct 2019 13.13920 12.96456 13.31384 12.87212 13.40628
## Nov 2019 13.11270 12.93661 13.28879 12.84340 13.38200
## Dec 2019 13.24990 13.07235 13.42745 12.97836 13.52144
## Jan 2020 13.12904 12.94947 13.30862 12.85440 13.40368
## Feb 2020 13.10545 12.92428 13.28661 12.82838 13.38251
## Mar 2020 13.19776 13.01510 13.38041 12.91841 13.47710
## Apr 2020 13.15239 12.96821 13.33657 12.87071 13.43407
## May 2020 13.12928 12.94361 13.31495 12.84532 13.41323
## Jun 2020 13.08366 12.89650 13.27081 12.79742 13.36989
## Jul 2020 13.16069 12.97207 13.34932 12.87221 13.44918
## Aug 2020 13.16595 12.97586 13.35604 12.87523 13.45667
## Sep 2020 13.11209 12.92055 13.30363 12.81915 13.40502
## Oct 2020 13.16245 12.96946 13.35543 12.86731 13.45759
## Nov 2020 13.13589 12.94148 13.33029 12.83857 13.43320
## Dec 2020 13.27290 13.07706 13.46873 12.97339 13.57241
## Jan 2021 13.15220 12.95439 13.35001 12.84968 13.45472
## Feb 2021 13.12849 12.92914 13.32784 12.82361 13.43337
## Mar 2021 13.22087 13.02006 13.42168 12.91376 13.52798
## Apr 2021 13.17566 12.97337 13.37796 12.86628 13.48504
## May 2021 13.15255 12.94880 13.35629 12.84094 13.46415
Como se puede ver, el modelo es bueno para predecir tendencias. Sin embargo, no se ve muy acertado para los datos que se tiene en el conjunto de test.
###Prophet algoritmo de FB
## Warning: package 'Rcpp' was built under R version 4.0.5## Warning: package 'rlang' was built under R version 4.0.5## Warning: package 'prophet' was built under R version 4.0.5## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Veamos que tanto se acerca la predicción a los datos reales, en el siguiente gráfico. Los datos del conjunto de prueba se representan en rojo mientras los datos de la predicción se representan en azul.
como podemos observar, las predicciones realizadas por prophet fueron peores y erroneas en comparacion con la funcion de Arima.
##Prediccion de los ultimos 3 años
train <- head(consumoGS_ts, round(length(consumoGS_ts) * 0.891632653))
h <- length(consumoGS_ts) - length(train)
test <- tail(consumoGS_ts, h)
end(train)## [1] 2019 1fitArima %>%
forecast(h) ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 13.03072 12.96888 13.09256 12.93614 13.12530
## Feb 2015 13.00463 12.93677 13.07249 12.90085 13.10841
## Mar 2015 13.09947 13.02944 13.16950 12.99237 13.20658
## Apr 2015 13.04365 12.96995 13.11735 12.93094 13.15636
## May 2015 13.02189 12.94554 13.09824 12.90512 13.13866
## Jun 2015 12.97006 12.89072 13.04940 12.84872 13.09139
## Jul 2015 13.04983 12.96784 13.13183 12.92443 13.17523
## Aug 2015 13.05805 12.97337 13.14273 12.92854 13.18755
## Sep 2015 13.01030 12.92308 13.09752 12.87690 13.14370
## Oct 2015 13.04912 12.95939 13.13885 12.91189 13.18635
## Nov 2015 13.02783 12.93569 13.11997 12.88691 13.16875
## Dec 2015 13.16829 13.07378 13.26281 13.02374 13.31284
## Jan 2016 13.04320 12.94569 13.14071 12.89407 13.19233
## Feb 2016 13.02564 12.92563 13.12564 12.87270 13.17858
## Mar 2016 13.11409 13.01176 13.21642 12.95758 13.27059
## Apr 2016 13.06203 12.95736 13.16669 12.90196 13.22209
## May 2016 13.03888 12.93197 13.14579 12.87537 13.20239
## Jun 2016 13.00052 12.89139 13.10965 12.83362 13.16743
## Jul 2016 13.07097 12.95968 13.18227 12.90076 13.24119
## Aug 2016 13.07571 12.96229 13.18914 12.90224 13.24919
## Sep 2016 13.02441 12.90889 13.13992 12.84774 13.20107
## Oct 2016 13.06679 12.94922 13.18437 12.88698 13.24661
## Nov 2016 13.04360 12.92401 13.16318 12.86071 13.22648
## Dec 2016 13.17061 13.04903 13.29219 12.98468 13.35655
## Jan 2017 13.05617 12.93330 13.17904 12.86826 13.24408
## Feb 2017 13.03266 12.90822 13.15709 12.84235 13.22296
## Mar 2017 13.12479 12.99868 13.25089 12.93193 13.31764
## Apr 2017 13.08114 12.95346 13.20883 12.88586 13.27642plot(forecast(fitArima,h))
df<-data.frame(ds=as.Date(as.yearmon(time(train))),y=as.matrix(train) )
testdf<-data.frame(ds=as.Date(as.yearmon(time(test))),y=as.matrix(test) )
fitProphet<-prophet(df,yearly.seasonality = T,weekly.seasonality = T)## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.future <- make_future_dataframe(fitProphet,periods = h,freq = "month", include_history = T)
p <- predict(fitProphet,future)
p<-p[,c("ds","yhat","yhat_lower","yhat_upper")]
plot(fitProphet,p)
pred<-tail(p,h)
pred$y<-testdf$y
ggplot(pred, aes(x=ds, y=yhat)) +
geom_line(size=1, alpha=0.8) +
geom_ribbon(aes(ymin=yhat_lower, ymax=yhat_upper), fill="blue", alpha=0.2) +
geom_line(data=pred, aes(x=ds, y=y),color="red")
En la grafica de los ultimos tres anios podemos ver que fue accurate su predccion con un bajon de consumo de gasolina a mediados de 2020. La prediccion dio lugar a que la tendencia fuera media y nomral, pero en los datos normal (en rojo) podemos ver que se tuvo un decrecimeinto significativo cuando se tuvo la pandemia, lo cual tiene sentido ya que las restricciones no permitian la movilizacion por carro a inicios de esta.
Serie de tiempo Gasolina Super Importacion
Union de columnas de tiempo
Creando la serie de tiempo. Viendo incio, fin y frencuencia
## [1] 2001 1## [1] 2021 5## [1] 12Creando el grafico de la serie
Se observa que es una serie de tiempo no estacionaria ya que su componente tendencia no es constante, si no que es creciente. Ademas, la varianza tampoco es constante.
Descomponiendo la serie
Se puede observar que la componente tendencia es creciente, la componente estacional si tiene un patron y la varianza no es constante. Se tendra que transformar
Transformacion logaritmica para estacionar serie
logSS <- log(trainS)
lambda <- BoxCox.lambda(trainS)
boxcoxGR<-BoxCox(trainS,lambda)
plot(decompose(logSS))
plot(logSS)
Al hacer la transformacion logaritmica, se observa que la serie se ha estabilizado un poco mas en media. Sin embargo en varianza sigue lejos de estar constante
Ahora se haran las pruebas de DickeyFuller
adfTest(trainS)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -0.567
## P VALUE:
## 0.4339
##
## Description:
## Fri Aug 06 21:43:01 2021 by user: jfdelunitrootTest(trainS)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: -0.567
## P VALUE:
## t: 0.4704
## n: 0.5542
##
## Description:
## Fri Aug 06 21:43:01 2021 by user: jfdelCon las pruebas de Dickey-Fuller, se observa que p es mayor a 0.05. No se puede rechazar Ho y existen raices unitarias. Aunque la serie se vea mas estable en media, aun no es estacionaria en media
Ahora se hace con una diferenciacion en la serie de tiempo
GS_diff<-diff(logSS)
adfTest(diff(trainS))## Warning in adfTest(diff(trainS)): p-value smaller than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -16.9539
## P VALUE:
## 0.01
##
## Description:
## Fri Aug 06 21:43:01 2021 by user: jfdelAhora el valor p si es menor a 0.05. Ahora si se puede rechazar la Ho, y ahora no hay raices unitarias
Valores de p y q
Se observa que el ACF se anula muy rapido y cerca de uno por lo que p puede ser uno. En el PACF, se observa que q es 1 ya que hay un pico en el retardo que se sale en 1. Por lo tanto tenemos un modelo ARIMA(1,1,1)
Valores p y q del la componente estacional
p es 1 ya que se tiene solamente un ciclo antes de que se anule q es 1 ya que se acerca a la frecuencia de 12, antes de anularse Por lo tanto, tenemos un modelo ARIMA (1,1,1)
Ahora se comprobara mediante fitarima si en efecto, los coeficiente son significativos
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.418542 0.077040 -5.4328 5.549e-08 ***
## ma1 -0.856811 0.045973 -18.6372 < 2.2e-16 ***
## sar1 -0.197184 0.080253 -2.4570 0.01401 *
## sma1 -0.999990 0.099506 -10.0495 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.654059 0.148938 -4.3915 1.126e-05 ***
## ma1 -0.537645 0.181707 -2.9589 0.003088 **
## ma2 -0.299850 0.177230 -1.6919 0.090670 .
## sma1 -0.161559 0.081102 -1.9921 0.046365 *
## drift 1299.928853 616.505243 2.1085 0.034984 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como se puede observar, los coeficientes son significativos.
Ahora se hace el analisis de los residuales
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,1,1)[12]
## Q* = 19.747, df = 20, p-value = 0.4738
##
## Model df: 4. Total lags used: 24Los residuos se comportan de manera normal y continua. Ademas, en la prueba L-jung-Box se tuvo un valor de p mayor a 0.05. Por lo que es un modelo apto para predecir.
Ahora se hace la prediccion arima en si.
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2015 13.15250 12.85453 13.45046 12.69680 13.60819
## Jun 2015 13.19007 12.88110 13.49905 12.71754 13.66261
## Jul 2015 13.05966 12.74122 13.37810 12.57264 13.54667
## Aug 2015 13.11468 12.79608 13.43328 12.62742 13.60194
## Sep 2015 13.04087 12.71997 13.36177 12.55010 13.53165
## Oct 2015 13.20716 12.88517 13.52914 12.71472 13.69959
## Nov 2015 13.10977 12.78623 13.43331 12.61496 13.60458
## Dec 2015 13.23458 12.90986 13.55929 12.73797 13.73118
## Jan 2016 13.17263 12.84683 13.49842 12.67437 13.67089
## Feb 2016 13.15992 12.83276 13.48709 12.65957 13.66028
## Mar 2016 13.23094 12.90240 13.55948 12.72848 13.73340
## Apr 2016 13.12845 12.79855 13.45836 12.62391 13.63300
## May 2016 13.24550 12.91553 13.57546 12.74086 13.75013
## Jun 2016 13.17113 12.83867 13.50360 12.66267 13.67960
## Jul 2016 13.14640 12.81334 13.47947 12.63702 13.65578
## Aug 2016 13.11720 12.78290 13.45149 12.60594 13.62846
## Sep 2016 13.12343 12.78819 13.45866 12.61073 13.63613
## Oct 2016 13.18349 12.84722 13.51976 12.66921 13.69777
## Nov 2016 13.17294 12.83565 13.51024 12.65709 13.68879
## Dec 2016 13.29044 12.95228 13.62861 12.77326 13.80762
## Jan 2017 13.19911 12.86005 13.53818 12.68056 13.71766
## Feb 2017 13.25628 12.91622 13.59633 12.73621 13.77635
## Mar 2017 13.30061 12.95957 13.64165 12.77904 13.82219
## Apr 2017 13.21080 12.86878 13.55283 12.68772 13.73388
## May 2017 13.27490 12.92861 13.62119 12.74530 13.80451
## Jun 2017 13.22261 12.87586 13.56936 12.69230 13.75292
## Jul 2017 13.17704 12.82870 13.52538 12.64430 13.70978
## Aug 2017 13.16445 12.81507 13.51382 12.63013 13.69877
## Sep 2017 13.15489 12.80426 13.50552 12.61865 13.69113
## Oct 2017 13.23590 12.88414 13.58766 12.69793 13.77386
## Nov 2017 13.20823 12.85524 13.56121 12.66839 13.74807
## Dec 2017 13.32717 12.97320 13.68114 12.78582 13.86852
## Jan 2018 13.24163 12.88654 13.59673 12.69856 13.78471
## Feb 2018 13.28502 12.92878 13.64126 12.74020 13.82984
## Mar 2018 13.33462 12.97724 13.69200 12.78806 13.88118
## Apr 2018 13.24231 12.88380 13.60082 12.69401 13.79060
## May 2018 13.31685 12.95568 13.67801 12.76449 13.86920
## Jun 2018 13.26020 12.89823 13.62217 12.70662 13.81379
## Jul 2018 13.21874 12.85532 13.58216 12.66294 13.77455
## Aug 2018 13.20287 12.83830 13.56745 12.64531 13.76044
## Sep 2018 13.19643 12.83058 13.56228 12.63691 13.75595
## Oct 2018 13.27331 12.90627 13.64034 12.71197 13.83464
## Nov 2018 13.24901 12.88069 13.61734 12.68571 13.81231
## Dec 2018 13.36767 12.99834 13.73701 12.80282 13.93252
## Jan 2019 13.28099 12.91040 13.65159 12.71421 13.84778
## Feb 2019 13.32710 12.95532 13.69888 12.75851 13.89569
## Mar 2019 13.37566 13.00269 13.74862 12.80526 13.94605
## Apr 2019 13.28384 12.90970 13.65798 12.71164 13.85604
## May 2019 13.35632 12.97916 13.73348 12.77950 13.93314
## Jun 2019 13.30054 12.92258 13.67849 12.72251 13.87857
## Jul 2019 13.25826 12.87875 13.63777 12.67785 13.83868
## Aug 2019 13.24304 12.86233 13.62375 12.66080 13.82528
## Sep 2019 13.23598 12.85392 13.61804 12.65167 13.82030
## Oct 2019 13.31367 12.93037 13.69698 12.72746 13.89989
## Nov 2019 13.28872 12.90405 13.67338 12.70042 13.87701
## Dec 2019 13.40743 13.02171 13.79314 12.81753 13.99733
## Jan 2020 13.32098 12.93386 13.70809 12.72893 13.91302
## Feb 2020 13.36655 12.97819 13.75490 12.77261 13.96048
## Mar 2020 13.41531 13.02572 13.80490 12.81948 14.01114
## Apr 2020 13.32339 12.93257 13.71422 12.72568 13.92111
## May 2020 13.39628 13.00244 13.79013 12.79395 13.99861
## Jun 2020 13.34033 12.94563 13.73502 12.73670 13.94396
## Jul 2020 13.29821 12.90192 13.69451 12.69213 13.90430
## Aug 2020 13.28286 12.88532 13.68041 12.67487 13.89086
## Sep 2020 13.27593 12.87697 13.67488 12.66578 13.88608
## Oct 2020 13.35346 12.95321 13.75370 12.74134 13.96558
## Nov 2020 13.32863 12.92696 13.73030 12.71433 13.94293
## Dec 2020 13.44733 13.04458 13.85008 12.83138 14.06329
## Jan 2021 13.36084 12.95657 13.76511 12.74256 13.97911
## Feb 2021 13.40651 13.00096 13.81207 12.78627 14.02675
## Mar 2021 13.45523 13.04840 13.86207 12.83303 14.07744
## Apr 2021 13.36334 12.95522 13.77145 12.73918 13.98749
## May 2021 13.43615 13.02495 13.84734 12.80728 14.06501
Se puede observar que el algoritmo pudo predecir la tendencia, sin embargo, no pudo predecir los datos.
Analisis con profet de facebook gasolina super importacion
## function (package, help, pos = 2, lib.loc = NULL, character.only = FALSE,
## logical.return = FALSE, warn.conflicts, quietly = FALSE,
## verbose = getOption("verbose"), mask.ok, exclude, include.only,
## attach.required = missing(include.only))
## {
## conf.ctrl <- getOption("conflicts.policy")
## if (is.character(conf.ctrl))
## conf.ctrl <- switch(conf.ctrl, strict = list(error = TRUE,
## warn = FALSE), depends.ok = list(error = TRUE, generics.ok = TRUE,
## can.mask = c("base", "methods", "utils", "grDevices",
## "graphics", "stats"), depends.ok = TRUE), warning(gettextf("unknown conflict policy: %s",
## sQuote(conf.ctrl)), call. = FALSE, domain = NA))
## if (!is.list(conf.ctrl))
## conf.ctrl <- NULL
## stopOnConflict <- isTRUE(conf.ctrl$error)
## if (missing(warn.conflicts))
## warn.conflicts <- if (isFALSE(conf.ctrl$warn))
## FALSE
## else TRUE
## if ((!missing(include.only)) && (!missing(exclude)))
## stop(gettext("only one of 'include.only' and 'exclude' can be used"),
## call. = FALSE, domain = NA)
## testRversion <- function(pkgInfo, pkgname, pkgpath) {
## if (is.null(built <- pkgInfo$Built))
## stop(gettextf("package %s has not been installed properly\n",
## sQuote(pkgname)), call. = FALSE, domain = NA)
## R_version_built_under <- as.numeric_version(built$R)
## if (R_version_built_under < "3.0.0")
## stop(gettextf("package %s was built before R 3.0.0: please re-install it",
## sQuote(pkgname)), call. = FALSE, domain = NA)
## current <- getRversion()
## if (length(Rdeps <- pkgInfo$Rdepends2)) {
## for (dep in Rdeps) if (length(dep) > 1L) {
## target <- dep$version
## res <- do.call(dep$op, if (is.character(target))
## list(as.numeric(R.version[["svn rev"]]), as.numeric(sub("^r",
## "", target)))
## else list(current, as.numeric_version(target)))
## if (!res)
## stop(gettextf("This is R %s, package %s needs %s %s",
## current, sQuote(pkgname), dep$op, target),
## call. = FALSE, domain = NA)
## }
## }
## if (R_version_built_under > current)
## warning(gettextf("package %s was built under R version %s",
## sQuote(pkgname), as.character(built$R)), call. = FALSE,
## domain = NA)
## platform <- built$Platform
## r_arch <- .Platform$r_arch
## if (.Platform$OS.type == "unix") {
## }
## else {
## if (nzchar(platform) && !grepl("mingw", platform))
## stop(gettextf("package %s was built for %s",
## sQuote(pkgname), platform), call. = FALSE,
## domain = NA)
## }
## if (nzchar(r_arch) && file.exists(file.path(pkgpath,
## "libs")) && !file.exists(file.path(pkgpath, "libs",
## r_arch)))
## stop(gettextf("package %s is not installed for 'arch = %s'",
## sQuote(pkgname), r_arch), call. = FALSE, domain = NA)
## }
## checkNoGenerics <- function(env, pkg) {
## nenv <- env
## ns <- .getNamespace(as.name(pkg))
## if (!is.null(ns))
## nenv <- asNamespace(ns)
## if (exists(".noGenerics", envir = nenv, inherits = FALSE))
## TRUE
## else {
## !any(startsWith(names(env), ".__T"))
## }
## }
## checkConflicts <- function(package, pkgname, pkgpath, nogenerics,
## env) {
## dont.mind <- c("last.dump", "last.warning", ".Last.value",
## ".Random.seed", ".Last.lib", ".onDetach", ".packageName",
## ".noGenerics", ".required", ".no_S3_generics", ".Depends",
## ".requireCachedGenerics")
## sp <- search()
## lib.pos <- which(sp == pkgname)
## ob <- names(as.environment(lib.pos))
## if (!nogenerics) {
## these <- ob[startsWith(ob, ".__T__")]
## gen <- gsub(".__T__(.*):([^:]+)", "\\1", these)
## from <- gsub(".__T__(.*):([^:]+)", "\\2", these)
## gen <- gen[from != package]
## ob <- ob[!(ob %in% gen)]
## }
## ipos <- seq_along(sp)[-c(lib.pos, match(c("Autoloads",
## "CheckExEnv"), sp, 0L))]
## cpos <- NULL
## conflicts <- vector("list", 0)
## for (i in ipos) {
## obj.same <- match(names(as.environment(i)), ob, nomatch = 0L)
## if (any(obj.same > 0L)) {
## same <- ob[obj.same]
## same <- same[!(same %in% dont.mind)]
## Classobjs <- which(startsWith(same, ".__"))
## if (length(Classobjs))
## same <- same[-Classobjs]
## same.isFn <- function(where) vapply(same, exists,
## NA, where = where, mode = "function", inherits = FALSE)
## same <- same[same.isFn(i) == same.isFn(lib.pos)]
## not.Ident <- function(ch, TRAFO = identity, ...) vapply(ch,
## function(.) !identical(TRAFO(get(., i)), TRAFO(get(.,
## lib.pos)), ...), NA)
## if (length(same))
## same <- same[not.Ident(same)]
## if (length(same) && identical(sp[i], "package:base"))
## same <- same[not.Ident(same, ignore.environment = TRUE)]
## if (length(same)) {
## conflicts[[sp[i]]] <- same
## cpos[sp[i]] <- i
## }
## }
## }
## if (length(conflicts)) {
## if (stopOnConflict) {
## emsg <- ""
## pkg <- names(conflicts)
## notOK <- vector("list", 0)
## for (i in seq_along(conflicts)) {
## pkgname <- sub("^package:", "", pkg[i])
## if (pkgname %in% canMaskEnv$canMask)
## next
## same <- conflicts[[i]]
## if (is.list(mask.ok))
## myMaskOK <- mask.ok[[pkgname]]
## else myMaskOK <- mask.ok
## if (isTRUE(myMaskOK))
## same <- NULL
## else if (is.character(myMaskOK))
## same <- setdiff(same, myMaskOK)
## if (length(same)) {
## notOK[[pkg[i]]] <- same
## msg <- .maskedMsg(sort(same), pkg = sQuote(pkg[i]),
## by = cpos[i] < lib.pos)
## emsg <- paste(emsg, msg, sep = "\n")
## }
## }
## if (length(notOK)) {
## msg <- gettextf("Conflicts attaching package %s:\n%s",
## sQuote(package), emsg)
## stop(errorCondition(msg, package = package,
## conflicts = conflicts, class = "packageConflictError"))
## }
## }
## if (warn.conflicts) {
## packageStartupMessage(gettextf("\nAttaching package: %s\n",
## sQuote(package)), domain = NA)
## pkg <- names(conflicts)
## for (i in seq_along(conflicts)) {
## msg <- .maskedMsg(sort(conflicts[[i]]), pkg = sQuote(pkg[i]),
## by = cpos[i] < lib.pos)
## packageStartupMessage(msg, domain = NA)
## }
## }
## }
## }
## if (verbose && quietly)
## message("'verbose' and 'quietly' are both true; being verbose then ..")
## if (!missing(package)) {
## if (is.null(lib.loc))
## lib.loc <- .libPaths()
## lib.loc <- lib.loc[dir.exists(lib.loc)]
## if (!character.only)
## package <- as.character(substitute(package))
## if (length(package) != 1L)
## stop("'package' must be of length 1")
## if (is.na(package) || (package == ""))
## stop("invalid package name")
## pkgname <- paste0("package:", package)
## newpackage <- is.na(match(pkgname, search()))
## if (newpackage) {
## pkgpath <- find.package(package, lib.loc, quiet = TRUE,
## verbose = verbose)
## if (length(pkgpath) == 0L) {
## if (length(lib.loc) && !logical.return)
## stop(packageNotFoundError(package, lib.loc,
## sys.call()))
## txt <- if (length(lib.loc))
## gettextf("there is no package called %s", sQuote(package))
## else gettext("no library trees found in 'lib.loc'")
## if (logical.return) {
## warning(txt, domain = NA)
## return(FALSE)
## }
## else stop(txt, domain = NA)
## }
## which.lib.loc <- normalizePath(dirname(pkgpath),
## "/", TRUE)
## pfile <- system.file("Meta", "package.rds", package = package,
## lib.loc = which.lib.loc)
## if (!nzchar(pfile))
## stop(gettextf("%s is not a valid installed package",
## sQuote(package)), domain = NA)
## pkgInfo <- readRDS(pfile)
## testRversion(pkgInfo, package, pkgpath)
## if (is.character(pos)) {
## npos <- match(pos, search())
## if (is.na(npos)) {
## warning(gettextf("%s not found on search path, using pos = 2",
## sQuote(pos)), domain = NA)
## pos <- 2
## }
## else pos <- npos
## }
## deps <- unique(names(pkgInfo$Depends))
## depsOK <- isTRUE(conf.ctrl$depends.ok)
## if (depsOK) {
## canMaskEnv <- dynGet("__library_can_mask__",
## NULL)
## if (is.null(canMaskEnv)) {
## canMaskEnv <- new.env()
## canMaskEnv$canMask <- union("base", conf.ctrl$can.mask)
## "__library_can_mask__" <- canMaskEnv
## }
## canMaskEnv$canMask <- unique(c(package, deps,
## canMaskEnv$canMask))
## }
## else canMaskEnv <- NULL
## if (attach.required)
## .getRequiredPackages2(pkgInfo, quietly = quietly)
## cr <- conflictRules(package)
## if (missing(mask.ok))
## mask.ok <- cr$mask.ok
## if (missing(exclude))
## exclude <- cr$exclude
## if (packageHasNamespace(package, which.lib.loc)) {
## if (isNamespaceLoaded(package)) {
## newversion <- as.numeric_version(pkgInfo$DESCRIPTION["Version"])
## oldversion <- as.numeric_version(getNamespaceVersion(package))
## if (newversion != oldversion) {
## tryCatch(unloadNamespace(package), error = function(e) {
## P <- if (!is.null(cc <- conditionCall(e)))
## paste("Error in", deparse(cc)[1L], ": ")
## else "Error : "
## stop(gettextf("Package %s version %s cannot be unloaded:\n %s",
## sQuote(package), oldversion, paste0(P,
## conditionMessage(e), "\n")), domain = NA)
## })
## }
## }
## tt <- tryCatch({
## attr(package, "LibPath") <- which.lib.loc
## ns <- loadNamespace(package, lib.loc)
## env <- attachNamespace(ns, pos = pos, deps,
## exclude, include.only)
## }, error = function(e) {
## P <- if (!is.null(cc <- conditionCall(e)))
## paste(" in", deparse(cc)[1L])
## else ""
## msg <- gettextf("package or namespace load failed for %s%s:\n %s",
## sQuote(package), P, conditionMessage(e))
## if (logical.return)
## message(paste("Error:", msg), domain = NA)
## else stop(msg, call. = FALSE, domain = NA)
## })
## if (logical.return && is.null(tt))
## return(FALSE)
## attr(package, "LibPath") <- NULL
## {
## on.exit(detach(pos = pos))
## nogenerics <- !.isMethodsDispatchOn() || checkNoGenerics(env,
## package)
## if (isFALSE(conf.ctrl$generics.ok) || (stopOnConflict &&
## !isTRUE(conf.ctrl$generics.ok)))
## nogenerics <- TRUE
## if (stopOnConflict || (warn.conflicts && !exists(".conflicts.OK",
## envir = env, inherits = FALSE)))
## checkConflicts(package, pkgname, pkgpath,
## nogenerics, ns)
## on.exit()
## if (logical.return)
## return(TRUE)
## else return(invisible(.packages()))
## }
## }
## else stop(gettextf("package %s does not have a namespace and should be re-installed",
## sQuote(package)), domain = NA)
## }
## if (verbose && !newpackage)
## warning(gettextf("package %s already present in search()",
## sQuote(package)), domain = NA)
## }
## else if (!missing(help)) {
## if (!character.only)
## help <- as.character(substitute(help))
## pkgName <- help[1L]
## pkgPath <- find.package(pkgName, lib.loc, verbose = verbose)
## docFiles <- c(file.path(pkgPath, "Meta", "package.rds"),
## file.path(pkgPath, "INDEX"))
## if (file.exists(vignetteIndexRDS <- file.path(pkgPath,
## "Meta", "vignette.rds")))
## docFiles <- c(docFiles, vignetteIndexRDS)
## pkgInfo <- vector("list", 3L)
## readDocFile <- function(f) {
## if (basename(f) %in% "package.rds") {
## txt <- readRDS(f)$DESCRIPTION
## if ("Encoding" %in% names(txt)) {
## to <- if (Sys.getlocale("LC_CTYPE") == "C")
## "ASCII//TRANSLIT"
## else ""
## tmp <- try(iconv(txt, from = txt["Encoding"],
## to = to))
## if (!inherits(tmp, "try-error"))
## txt <- tmp
## else warning("'DESCRIPTION' has an 'Encoding' field and re-encoding is not possible",
## call. = FALSE)
## }
## nm <- paste0(names(txt), ":")
## formatDL(nm, txt, indent = max(nchar(nm, "w")) +
## 3L)
## }
## else if (basename(f) %in% "vignette.rds") {
## txt <- readRDS(f)
## if (is.data.frame(txt) && nrow(txt))
## cbind(basename(gsub("\\.[[:alpha:]]+$", "",
## txt$File)), paste(txt$Title, paste0(rep.int("(source",
## NROW(txt)), ifelse(nzchar(txt$PDF), ", pdf",
## ""), ")")))
## else NULL
## }
## else readLines(f)
## }
## for (i in which(file.exists(docFiles))) pkgInfo[[i]] <- readDocFile(docFiles[i])
## y <- list(name = pkgName, path = pkgPath, info = pkgInfo)
## class(y) <- "packageInfo"
## return(y)
## }
## else {
## if (is.null(lib.loc))
## lib.loc <- .libPaths()
## db <- matrix(character(), nrow = 0L, ncol = 3L)
## nopkgs <- character()
## for (lib in lib.loc) {
## a <- .packages(all.available = TRUE, lib.loc = lib)
## for (i in sort(a)) {
## file <- system.file("Meta", "package.rds", package = i,
## lib.loc = lib)
## title <- if (nzchar(file)) {
## txt <- readRDS(file)
## if (is.list(txt))
## txt <- txt$DESCRIPTION
## if ("Encoding" %in% names(txt)) {
## to <- if (Sys.getlocale("LC_CTYPE") == "C")
## "ASCII//TRANSLIT"
## else ""
## tmp <- try(iconv(txt, txt["Encoding"], to,
## "?"))
## if (!inherits(tmp, "try-error"))
## txt <- tmp
## else warning("'DESCRIPTION' has an 'Encoding' field and re-encoding is not possible",
## call. = FALSE)
## }
## txt["Title"]
## }
## else NA
## if (is.na(title))
## title <- " ** No title available ** "
## db <- rbind(db, cbind(i, lib, title))
## }
## if (length(a) == 0L)
## nopkgs <- c(nopkgs, lib)
## }
## dimnames(db) <- list(NULL, c("Package", "LibPath", "Title"))
## if (length(nopkgs) && !missing(lib.loc)) {
## pkglist <- paste(sQuote(nopkgs), collapse = ", ")
## msg <- sprintf(ngettext(length(nopkgs), "library %s contains no packages",
## "libraries %s contain no packages"), pkglist)
## warning(msg, domain = NA)
## }
## y <- list(header = NULL, results = db, footer = NULL)
## class(y) <- "libraryIQR"
## return(y)
## }
## if (logical.return)
## TRUE
## else invisible(.packages())
## }
## <bytecode: 0x000000001309d8c8>
## <environment: namespace:base>## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Como se puede observar, el ARIMA tuvo una mejor prediccion.
Prediccion ultimos 3 anios
## [1] 2019 2## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2015 13.15250 12.85453 13.45046 12.69680 13.60819
## Jun 2015 13.19007 12.88110 13.49905 12.71754 13.66261
## Jul 2015 13.05966 12.74122 13.37810 12.57264 13.54667
## Aug 2015 13.11468 12.79608 13.43328 12.62742 13.60194
## Sep 2015 13.04087 12.71997 13.36177 12.55010 13.53165
## Oct 2015 13.20716 12.88517 13.52914 12.71472 13.69959
## Nov 2015 13.10977 12.78623 13.43331 12.61496 13.60458
## Dec 2015 13.23458 12.90986 13.55929 12.73797 13.73118
## Jan 2016 13.17263 12.84683 13.49842 12.67437 13.67089
## Feb 2016 13.15992 12.83276 13.48709 12.65957 13.66028
## Mar 2016 13.23094 12.90240 13.55948 12.72848 13.73340
## Apr 2016 13.12845 12.79855 13.45836 12.62391 13.63300
## May 2016 13.24550 12.91553 13.57546 12.74086 13.75013
## Jun 2016 13.17113 12.83867 13.50360 12.66267 13.67960
## Jul 2016 13.14640 12.81334 13.47947 12.63702 13.65578
## Aug 2016 13.11720 12.78290 13.45149 12.60594 13.62846
## Sep 2016 13.12343 12.78819 13.45866 12.61073 13.63613
## Oct 2016 13.18349 12.84722 13.51976 12.66921 13.69777
## Nov 2016 13.17294 12.83565 13.51024 12.65709 13.68879
## Dec 2016 13.29044 12.95228 13.62861 12.77326 13.80762
## Jan 2017 13.19911 12.86005 13.53818 12.68056 13.71766
## Feb 2017 13.25628 12.91622 13.59633 12.73621 13.77635
## Mar 2017 13.30061 12.95957 13.64165 12.77904 13.82219
## Apr 2017 13.21080 12.86878 13.55283 12.68772 13.73388
## May 2017 13.27490 12.92861 13.62119 12.74530 13.80451
## Jun 2017 13.22261 12.87586 13.56936 12.69230 13.75292
## Jul 2017 13.17704 12.82870 13.52538 12.64430 13.70978
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Se puede observar que la prediccion para el 2021 no es aceptable. Ademas, por el efecto de la pandemia y su comportamiento irregular, la prediccion tuvo peores resultados despues de la misma.
####Linea de tiempo para consumo de gasolina Regular Importacion
consumoGR_ts <- xts(consumoRegular$cantRegular, consumoRegular$regularDate)
consumoGR_ts <-ts_ts(consumoGR_ts)
#Saber cuando empieza la serie y cuando termina
start(consumoGR_ts)## [1] 2000 1end(consumoGR_ts)## [1] 2021 5#Saber la frecuencia de la serie
frequency(consumoGR_ts)## [1] 12#Ver el gráfico de la serie
plot(consumoGR_ts)
###Analisis Importacion Gasolina Regular y Diesel
datosImport <- read.csv(file = 'DatosImportacionCombustibles.csv')
datosImport$Dia <- 1
datosImport$Date <- as.Date(with(datosImport, paste(Anio, Mes, Dia,sep="-")), "%Y-%m-%d")
cantRegular <- c(datosImport$GasolinaRegular)
regularDate <- c(datosImport$Date)
importRegular <- data.frame(cantRegular,regularDate)
cantDiesel <- c(datosImport$TotalDiesel)
dieselDate <- c(datosImport$Date)
importDiesel <- data.frame(cantDiesel,dieselDate)Inicio y fin de la importacion Gasolina Regular
importRegular_ts <- xts(importRegular$cantRegular, importRegular$regularDate)
importRegular_ts <-ts_ts(importRegular_ts)
start(importRegular_ts)## [1] 2001 1end(importRegular_ts) ## [1] 2021 5train <- head(importRegular_ts, round(length(importRegular_ts) * 0.7))
h <- length(importRegular_ts) - length(train)
test <- tail(importRegular_ts, h)Frecuencia serie importacion Gasolina Regular
frequency(importRegular_ts)## [1] 12Grafico de la serie importacion Gasolina Regular
plot(importRegular_ts)
abline(reg=lm(importRegular_ts~time(importRegular_ts)), col=c("red"))
De la grafica puedo observar que tiene tendencia a aumentar. No es estacionaria media
Descomposicion serie importacion Gasolina Regular
dec.GS<-decompose(importRegular_ts)
plot(dec.GS)
Por la descomposicion tampoco parece ser estacionaria en varianza. El factor aleatorio parece tener un patron por lo que no sera una serie util para predecir. Es necesario transformarla
Transformacion log y grafico autocorrelacion serie importacion Gasolina Regular
logIR <- log(train)
lambda <- BoxCox.lambda(train)
boxcoxGR<-BoxCox(train,lambda)
plot(decompose(logIR))
plot(logIR)
acf(logIR)
No hay estacionariedad media
Prueba Dickey-Fuller
adfTest(train)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -0.6179
## P VALUE:
## 0.4177
##
## Description:
## Fri Aug 06 21:43:06 2021 by user: jfdelunitrootTest(train)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: -0.6179
## P VALUE:
## t: 0.4484
## n: 0.5437
##
## Description:
## Fri Aug 06 21:43:06 2021 by user: jfdelEl test Dickey-Fuller resulta con valores mayores a 0.05 por lo que Ho no se rechaza. Si hay raices unitarias. Se debe hacer una diferenciacion de la serie de tiempo.
Diferenciacion serie importacion Gasolina Regular
IR_diff<-diff(logIR)
adfTest(diff(train))## Warning in adfTest(diff(train)): p-value smaller than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -15.4037
## P VALUE:
## 0.01
##
## Description:
## Fri Aug 06 21:43:06 2021 by user: jfdelAl hacer este cambio el valor p si es menor a 0.05, se rechaza Ho y ahora no hay raices unitarias
Valores p y q
acf(logIR,50) 
pacf(logIR,50) 
Inicio y fin de la importacion Gasolina Diesel
importDiesel_ts <- xts(importDiesel$cantDiesel, importDiesel$dieselDate)
importDiesel_ts <-ts_ts(importDiesel_ts)
start(importDiesel_ts)## [1] 2001 1end(importDiesel_ts) ## [1] 2021 5train <- head(importDiesel_ts, round(length(importDiesel_ts) * 0.7))
h <- length(importDiesel_ts) - length(train)
test <- tail(importDiesel_ts, h)Frecuencia serie importacion Gasolina Diesel
frequency(importDiesel_ts)## [1] 12Grafico de la serie importacion Gasolina Diesel
plot(importDiesel_ts)
abline(reg=lm(importDiesel_ts~time(importDiesel_ts)), col=c("red"))
De la grafica puedo observar que tiene tendencia a aumentar. No es estacionaria media
Descomposicion serie importacion Gasolina Diesel
dec.GD<-decompose(importDiesel_ts)
plot(dec.GD)
Por la descomposicion tampoco parece ser estacionaria en varianza. El factor aleatorio parece tener un patron por lo que no sera una serie util para predecir. Es necesario transformarla
Transformacion log y grafico autocorrelacion serie importacion Gasolina Diesel
logID <- log(train)
lambda <- BoxCox.lambda(train)
boxcoxGD<-BoxCox(train,lambda)
plot(decompose(logID))
plot(logID)
acf(logID)
No hay estacionariedad media
Prueba Dickey-Fuller
adfTest(train)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -0.9005
## P VALUE:
## 0.3276
##
## Description:
## Fri Aug 06 21:43:06 2021 by user: jfdelunitrootTest(train)##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: -0.9005
## P VALUE:
## t: 0.3248
## n: 0.4892
##
## Description:
## Fri Aug 06 21:43:06 2021 by user: jfdelEl test Dickey-Fuller resulta con valores mayores a 0.05 por lo que Ho se rechaza. No hay raices unitarias. Se debe hacer una diferenciacion de la serie de tiempo.
Diferenciacion serie importacion Gasolina Diesel
ID_diff<-diff(logID)
adfTest(diff(train))## Warning in adfTest(diff(train)): p-value smaller than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## Dickey-Fuller: -15.2594
## P VALUE:
## 0.01
##
## Description:
## Fri Aug 06 21:43:06 2021 by user: jfdelAl hacer este cambio el valor p si es menor a 0.05, ya se tienen raices unitarias.
Valores p y q
acf(logID,50) 
pacf(logID,50) 
Analisis Consumo
## Warning: package 'rstan' was built under R version 4.0.5## Loading required package: StanHeaders## Warning: package 'StanHeaders' was built under R version 4.0.5## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file##
## Attaching package: 'rstan'## The following object is masked from 'package:epiDisplay':
##
## lookupAnálisis exploratorio - Box-Plot Diesel, Regular
Luego de realizar el analisis exploratorio de las variables continuas, en relacion a los tipos de gasolina Diesel, Regular y Super, se observo que estas tienen una distribucion normal en base a los diagramas de caja mostrados anteriormente. En ciertos diagramas se puede apreciar sesgos positivos o negativos, sin embargo, estos no son representativos para indicar que no muestran una distribucion normal. En base a lo anterior, se puede indicar que se puede obtener una serie de tiempo del conjunto de datos. Es necesario determinar si la serie de tiempo del conjunto de datos presenta estacionalidad, asi como tambien determinar si se debe realizar una normalizacion en los datos para poder brindar una mejor prediccion que se lograra con los modelos.
Series de Tiempo
Al encontrar las series de tiempo se obtuvieron los datos de las graficas anteriores. Sin embargo es necesario descomponer estas:
Al descomponer la serie de tiempo se puede observar que ninguna de las series son estacionarias ya que la media no es constante al igual que la varianza.
Estacionariedad
Para deterimnar la estacionariedad se realizaron los siguientes graficos que muestra la tendencia de cada serie, factor aleatorio y estacionariedad.
En las graficas anteriores se puede apreicar como estas siguen un patron, esto es indicador de que cumplen con estacionariedad. Se podria indicar que la serie de tiempo del diesel es estacionaria en media, ya que se mantiene relativamente constante, sin embargo, en los ultimos años esta oresenta una baja. Mientras que la serie de tiempo de la gasolina regular no es estacionaria en media.
## Warning in adf.test(logDieselST): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: logDieselST
## Dickey-Fuller = -5.6671, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary##
## Augmented Dickey-Fuller Test
##
## data: logRegularST
## Dickey-Fuller = -1.8568, Lag order = 6, p-value = 0.6361
## alternative hypothesis: stationary
Modelos
Modelo Diesel
Modelo Regular
Prophet
Gasolina diesel
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Gasolina Regular
## ds y
## 1 2001-01-01 202645.2
## 2 2001-02-01 205531.0
## 3 2001-03-01 229499.6
## 4 2001-04-01 210680.4
## 5 2001-05-01 208164.3
## 6 2001-06-01 195088.7## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
En base a las graficas anteriores se puede indicar que el modelo fue eficiente Dado que se ajusto el conjunto de datos lo mas posible, se puede observar que hay una cercania entre la prediccion y los datos reales. El mejor modelo producido es el de diesel ya que tiene menos variacion, mientras que la gasolina regular se dificulto la prediccion.